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Abstract 

We present a new method to introduce phase-space fluctuations in transport the- 
ories, corresponding to a fuh implementation of the Boltzmann-Langevin equation 
for fermionic systems. It is based on the procedure originally developed by Bauer et 
al. for transport codes employing the test particle method. In the new procedure, 
the Pauli principle is carefully checked, leading to a good reproduction of the correct 
fluctuations in the "continuum limit" (/i — s- 0). Accurate tests are carried out in one 
and two dimensional idealized systems, and finally results for a full 3D application 
are shown. We stress the reliability of this method, which can be easily plugged into 
existing tranport codes using test particles, and its general applicability to systems 
characterized by instabilities, like for instance multifragmentation processes. 

Key words: Fluctuations; stochastic collision integral; fermionic systems; 
transport theories. 
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1 Introduction 



In recent years, the dynamics of heavy ion collisions at intermediate energy 
has been extensively investigated within the framework of transport theo- 
ries, such as the Nordheim approach, in which the Vlasov equation for the 
one-body phase space density, /(r, p, t), is supplemented with a Pauli-blocked 
Boltzmann collision term p!p] . The basic ingredients that enter the resulting 
transport equation, often called Boltzmann-Uehling-Uhlenbeck (BUU) equa- 
tion, are the self-consistent mean-field potential and the two-body scattering 
cross sections. These transport models hence describe the time evolution of the 
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reduced one-body density in phase-space and, consequently, they are suited 
for the description of one-body observables, such as inclusive particle spectra 
in nuclear collisions, average collective flows and excitations. However, they 
cannot provide a reliable description of fluctuation phenomena, such as multi- 
fragmentation processes, i.e. the break-up of excited nuclear systems into many 
pieces. In fact, neither fluctuations of one-body observables nor many-body 
correlations can be addressed with this class of mean-field models. Hence suit- 
able extensions, including fluctuations of the one-body density, have to be 
considered. 



An intense theoretical work on fluctuations in nuclear dynamics has started 
in the past years, also stimulated by the availability of large amounts of ex- 
perimental data on fragment formation in intermediate energy heavy ion colli- 
sions and the possibility to observe a liquid-gas phase transition piHI5]|5]|71ISll^ . 
In order to introduce fluctuations in transport theories, a number of dif- 
ferent avenues have been taken, that can be essentially reconducted to two 
different classes of models. One is the class of molecular dynamics models 
Pll0|ll|12|13|14j while the other kind is represented by stochastic mean-field 
approaches 



In molecular dynamics models the many-body state is represented by a simple 
product wave function, with or without antisymmetrization. The single par- 
ticle wave functions are assumed to have a fixed Gaussian shape. In this way, 
though nucleon wave functions are supposed to be independent (mean-field 
approximation), the use of localized wave packets induces many-body correla- 
tions both in mean-field propagation and hard two body scattering (collision 
integral), which is treated stochastically. Hence this way to introduce many- 
body correlations and produce a trajectory branching is essentially based on 
the use of empirical gaussian wave packets. If wave functions were allowed to 
assume any shape, the method would become identical to standard mean-field 
descriptions. While the wave function localization appears appropriate to de- 
scribe final fragmentation channels, where each single particle wave function 
should be localized within a fragment, the use of fixed shape localized wave 
packets in the full dynamics could affect the correct description of one-body 
effects, such as spinodal instabilities and zero sound propagation 

On the other side, in the so-called stochastic mean-field approaches, the stochas- 
tic extension of the transport treatment for the one-particle density is obtained 
by introducing a stochastic term representing the fluctuating part of the colli- 
sion integral [TBIT7] , in close analogy with the Langevin equation for a brown- 
ian motion. This can be derived as the next-order correction, in the equation 
describing the time evolution of /, with respect to the standard average col- 
lision integral, leading to the Boltzmann-Langevin (BL) equation. Thus, the 
system is still described solely in terms of the reduced one-body density /, 
but this function experiences a stochastic time evolution in response to the 
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random effect of the ffuctuating collision term. In this way density fluctuations 
are introduced, that are amplified when instabilities or bifurcations occur in 
the dynamics. This procedure is suitable also for addressing multifragmenta- 
tion phenomena, since fragments can be associated with the regions where the 
spacial density becomes larger, which finally can be reconstructed by sampling 
the one-body distribution function. 

A specific method for solving the Boltzmann-Langevin equation by direct nu- 
merical simulation was introduced in Refs. [TSfTU] . In this numerical implemen- 
tation, the one-particle density /(r, p) is represented on a lattice of grid points 
in phase space and the collision integral is treated by considering all possible 
transitions between phase space cells, adding a noise term whose features are 
related to the average rate of transitions between two specified initial cells 
and two final cells. The numerical implementation of this method has only 
been possible in two dimensions (2D) because it requires too large computer 
resources in 3D. 

Hence several approximate solutions of the BL equations have been formu- 
lated, mostly based on the projection of the BL noise only on a given dynam- 
ical variable (such as the local quadrupole tensor of the momentum distribu- 
tion) or on r space only [20f2T] . More tractable fluctuating terms, such as a 
stochastic force added to the mean-field potential, have also been proposed 
and extensively applied to multifragmentation studies [22][71l8] . 

However, the implementation of the full structure in phase space of the orig- 
inal BL term can still be considered as an important goal to reach. In fact, 
this would allow to treat a more general class of phenomena, where the cor- 
rect description of fluctuations and correlations in p space is essential (such 
as particle production and fragment velocity correlations for instance). More- 
over, also in the multi-fragmentation mechanism, that is dominated by spacial 
density fluctuations, a more accurate representation of the full phase space dy- 
namics, including fluctuations, would allow to improve the description of the 
fragment kinematical properties and correlations. 

We also stress the general interest of this effort. Indeed transport phenomena 
occur in many physical systems, for which a more precise description of the 
time evolution of the one-body distribution function, including the effect of 
many-body correlations, would be important. 

A flrst attempt to introduce a fluctuating collision term in a 3D transport 
approach was made by Bauer et al. [23] • This method can be implemented 
relatively easily into standard transport codes that adopt the scattering of 
pseudo-particles (or test particles) as a method of solution of the collision in- 
tegral and it consists essentially in forcing similar two-body collisions to occur 
for neighboring test particles, deflned according to a given distance in phase 
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space, so that effectively two nucleons are involved in each particular collision 
event. The distance should reproduce the phase-space shape of the nucleon 
wave function. In this way the random nature of the two-body scattering, that 
in the standard codes applies only to test particles and is washed-out when 
using a huge number of them, is transferred to entire nucleons. However, in 
the procedure proposed in Ref . [23j the Pauli blocking is checked only for the 
collision of the two original test particles and not for the entire swarm affected, 
leading to some unpleasant features in the description of fermionic systems. 
Indeed the Pauli-blocking violation introduces important inaccuracies in the 
fluctuations of the one-body density. 

In the present manuscript we present a new method to reconstruct the phase 
space nucleon wave function in mean-fleld approaches, in such a way that the 
Pauli-blocking is checked for the entire cloud of moved particles. This will im- 
prove the description of the fluctuation variance, approaching the one expected 
for fermionic systems. We will also pay special attention to the deflnition of 
the phase-space metric that would optimize the value of this variance. 

The paper is organized as follows: We will flrst recall the main ingredients 
of the BL theory, in order to connect the formalism with the numerical im- 
plementation adopted (Section 2). Then we will discuss in more detail the 
methods that have been proposed so far to solve the Boltzmann-Langevin 
equation (Section 3). The new procedure that we foHow to build fluctuations 
is presented in Section 4. Several results demonstrating and analyzing in detail 
the method are discussed from Section 5 to 8. Conclusions and perspectives 
are drawn in Section 9. 



2 The Boltzmann-Langevin equation 

Within the semi-classical framework, the stochastic transport equation of mo- 
tion for the one-body distribution function / can be expressed in the following 
form, 

/ = |/ - {h[f]J} = K[f] = K[f] + 5K[f] , (1) 

where the left side describes the coUisionless propagation of the individual 
particles in their common self-consistent one-body fleld, while the right side 
expresses the effect of the residual binary collisions. 

The collision term K[f] has a stochastic character. For example, the distance 
a particle travels in the medium before colliding is stochastic, as is the re- 
sulting scattering angle. At the Boltzmann level, Eq.([T]) includes only the 
average part of the collision term, K[f]. Usually, in nuclear systems, the quan- 
tum statistics is taken into account by adding suitable Fermi blocking in the 
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single-particle final states, leading to the Uehling-Uhlenbeck collision term 
pips] . The resulting nuclear Boltzmann equation exists in many implemen- 
tations that differ with respect to both the physical input (such as the types 
of constituents included, the form of their effective Hamiltonian, and their 
differential interaction cross sections) and the numerical methods employed 
(whether of test particle or lattice type) and various names have been employed 
in the literature, including BUU (for Boltzmann- Uehling-Uhlenbeck), VUU 
(for Vlasov-Uehling-Uhlenbeck), Landau- Vlasov and BNV (for Boltzmann- 
Nordheim-Vlasov) p!P|26|27ll2S] . The Boltzmann-Langevin treatment includes 
also the stochastic part of the collision term, ^i^[/] [TBl[T7] . 

In the simple physical scenario where the residual interaction can be considered 
as binary collisions that are well localized in space and time, the average part 
is given by [2HI25] : 

Kir, pi) = W{12; 34) [/1/2/3/4 - /1/2/3/4] , (2) 

234 

where fi is a short-hand notation for /(r, pj, t) and / = 1 — / is the associated 
Fermi blocking factor, g is the degeneracy factor. 

The basic transition rate is simply related to the differential cross section for 
the corresponding two-body scattering process, 

(12; 34) = v,2 ( ^] S{p, + P2 - Ps - P4) , (3) 



12^34 



being V12 = |vi — V2I the relative velocity, and it thus has corresponding 
symmetry properties, 

iy(12;34) = W{21;3A) = 1^(34; 12) . (4) 



Since it arises from the same elementary two-body processes, the stochastic 
part of the collision term is fully determined by the basic transition rate as 
well, as a manifestation of the fluctuation-dissipation theorem. With the col- 
lisions assumed to be local in space and time, the correlation function for the 
fluctuating part of the collision term has the following form, 

^ SK{r, pi, t) 6K{r', pi', t') y= C(pi, p/, r, t) 6{r - r') 6{t - t') , (5) 

where -< ■ >- denotes the average with respect to the ensemble of possible 
trajectories resulting from the current one-body density /. Furthermore, for 
elastic scattering, the correlation kernel is given by |16j : 

C(pi, pi', r, t) = 6iv E W^(12; 34)F(12; 34) 

234 
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+ ^[Vr(ll';34)F(ll';34) -2Vr(13;l'4)F(13;l'4)] , (6) 

34 



with the short-hand notations 6u' = h^6{pi — p[) and F{12; 34) = f 1/2/3/4 + 
/1/2/3/4 {D is the dimension of the space considered). The symmetry proper- 
ties (HI) of the transition rate ensure that the following sum rules hold, 



^C(pi,p2,r,t) = ^C(pi,p2,r,t) =0, (7) 
1 2 

J2 C{pi, P2, r, t)pi = <^(Pi' P2' t)p2 = , (8) 
1 2 

Zl'^(Pi'P2' J''''^)^! = Zl'^(Pi'P2'r,t)e2 = 0, (9) 

1 2 

where = pj/2m is the kinetic energy for a specified momentum. These sum 
rules express the fact that each of the elementary binary collisions conserves 
particle number, momentum, and energy, respectively. 



3 Methods to solve the BL equation 



3.1 Lattice calculations 



The numerical implementation of the Boltzmann-Langevin equation is accom- 
plished by correctly simulating the basic stochastic process, i.e. the stochastic 
transition rate among phase-space cells. One can realize that this point is 
very delicate in finite systems, when we work with a relatively small number 
of nucleons, but still need to build a smooth distribution function. There are 
different ways to overcome this difficulty. 

As mentioned in the Introduction, one possibility, proposed in Ref.[IB], is to 
solve the BL equation on a lattice; phase space is therefore divided into a 
number of cells, each one having volume As = ArAp (s denotes a point in 
phase space: s = (r,p)). 

Each transition involves four locations, and the collision integral arises from 
the sum of the transitions evaluated for all possible combinations of the cells: 
from this consideration one recognizes that the practicality of the method is 
limited by the huge computational effort required. In fact, only two-dimensional 
implementations exist [T81IT9] . Each transition represents a basic stochastic 
process, and, following the BL theory, the actual number of such transitions 
in a time step At is dispersed around the mean value according to a Poisson 
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distribution, so its variance amounts to: 




(10) 



where the average is given by: 



z>(12;34) 



Asi As2 Ass As4 



/i/2/3/4W^(12;34) 






■ 5(ri - r2)(5(ri 



r3)(5(r2 - T4)At. 



(11) 



This is a compact form of the fluctuation-dissipation theorem, since it predicts 
that the fluctuations are simply related to the mean number of dissipative pro- 
cesses. Then, a noise term 5z/(12; 34) is added to the mean number u{12] 34), 
and thus the actual number of transitions is a random number picked from a 
normal distribution, the center of which is the mean value, and where the width 
is given by Eq. ffTOl) . In refs. [T9||29] this procedure is proven to yield correct 
results for a 2D system of fermions interacting through hard two-body scatter- 
ing. Starting from a non-equilibrium situation (two touching Fermi spheres) 
the method is successful in reproducing the expected fluctuations, preserv- 
ing also the average trajectory of the system. In particular, the fluctuations 
introduced build the expected statistical value for the equilibrium one-body 
density variance: 



evaluated considering phase-space cells of volume As = h , for which the 
fluctuating transition rate is implemented. Here feq{s) denotes the equilibrium 

value of the one-body distribution function and Sf{s) = f{s) < /(s) >~. Of 

course, Eq. fll2l) holds also for volumes larger than and fluctuations are 
scaled accordingly. Also the co-variance, i.e. the correlation between density 
fluctuations in two different phase-space points, s and s', is well reproduced 



The same authors of Ref. |19j also verified that, for the same idealized system, 
but prepared at low spatial density, early fluctuations developed in momentum 
space are subsequently transmitted into density fluctuations and amplifled 
by the nuclear mean fleld [30], leading to large instabilities and a statistical 
population of fragments. 

3.2 The pseudo-particle correlation method 

The method of Ref. [18] was originally developed as a way to overcome the 
problems arising from the solution of transport equations with the test parti- 
cle method, which we now discuss briefly. Within the test particle method, each 



a%s) =^ 6fis)5fis) y= Us)Us)h^/{gAs) 



(12) 
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nucleoli is represented by a collection of Ntest test particles, that are propa- 
gated according to the mean-field interaction and random two-body scattering. 
The use of a large number of test particles allows to have a smooth distribution 
function and a good coverage of the phase space for the Pauli blocking. On the 
other side, since collisions are treated stochastically for the single test parti- 
cles, the dispersion around the average number of nucleon collisions (Eq. fllOl) ) 
is automatically divided by Ntest- Hence fiuctuations which are introduced by 
the test particle algorithm can also be seen as the correct ones expected from 
the BL approach, divided by the factor Ntest- 

Various paths have been followed to overcome this problem. For instance, one 
can choose Ntest = 1? i-e. work with whole nucleons [TU], but in this case one 
has to solve the problem of the induced numerical errors on the smooth path 
of the dynamics, i.e. the mean-field and the Pauli-blo eking factors. Alterna- 
tively, Bauer et al. proposed a method to introduce a correlation between close 
particles in phase-space [25] . 

The method follows the idea, applied in extended Time-Dependent-Hartree- 
Fock (TDHF) calculations, of evolving the single-particle density including a 
statistical mixing of Slater determinants [31]. The jump from one configuration 
to the other is possible under the important assumption that any coherence be- 
tween determinants can be ignored (decorrelation approximation). The choice 
of the single-particle basis for the determinant is however somewhat arbitrary, 
apart from the requirement of momentum-energy conservation. These ideas 
can be translated into the semi-classical approximation using the test particle 
method. The mixing between Slater determinants is realized in the collision 
integral by means of the following procedure. 

(1) First of all, the nucleon- nucleon (NN) cross section is reduced by a factor 

Ntest- 

(2) Then two test particles ii and 12 are chosen as colliding partners, and will 
be moved from their positions pi and p2, to new positions pa and p4, 
according to the corresponding transition probability including the Pauli 
blocking of the final states [2]. 

(3) If the two particles can collide, the scattering actually involves two "clouds" 
of neighbouring test particles, corresponding to two entire nucleons (2 x 
Ntest particles). The contiguity criterion is based on the following defini- 
tion of phase-space distance: 

<■ = (p. - P.) VpI + (r. - r,) Vi?' (13) 

where p_F is the Fermi momentum and R is the radius of the considered 
system. (It should be noticed that the choice of this phase-space metric 
is rather arbitrary). 

(4) In order to ensure energy and momentum conservation, the final states 
are adjusted using the average momenta of the two clouds as initial states. 
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Two important remarks have to be done: 



• the Pauh blocking is checked only for the collision of the test particles ii 
and ^2 and not for each particle of the cloud. We will see in the following 
that this choice induces a strong violation of the Pauli principle; 

• the effect of the collision is a mere translation (in momentum space) of the 
two clouds to final positions, as sketched in Fig. [TJ 



Fig. 1. Illustrative picture of nucleon-nucleon scattering in the pseudo-particle cor- 
relation method. 

This method is conceptually simple, and moreover the implementation into 
existing BUU or BNV codes is really straightforward. Indeed it was immedi- 
ately used to test the effects of the fluctuations on fragment mass spectra in 
a 20iVe + 20iVe reaction at WOAMeV [23|. 

Its validity has been checked later by Chapelle et al. [29], for the same ideal- 
ized 2D system for which the BL lattice simulation method was implemented. 
Accordingly, a different metric, involving only distances in momentum space, 
was used. As shown in |29], the final equilibrium value obtained for the fluc- 
tuation variance does not reproduce the expected one. The variance profile, 
as a function of the energy, appears to be proportional to f^q rather than to 
feqfeq- Actually, this would be the result for a classical system, obeying to 
the Boltzmann statistics. This can be considered as a consequence of the fact 
that Pauli blocking is checked only for the first two colliding partners, and 
the system progressively evolves towards an equilibrium configuration which 
is more consistent with the one of a classical system. 

From this discussion it is evident that a proper treatment of the Pauli block- 
ing is a fundamental requirement for any numerical implementation of the 
Boltzmann-Langevin theory for fermionic systems. However, apart from these 
problems, the pseudo-particle correlation method represents a simple and 
practical way to implement fluctuations into transport codes. Therefore, it 
is worthwhile to think about some improvements to make it more accurate. 
This is the aim of the present work, as we will discuss in the next Section. 
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4 The improved pseudo-particle correlation method 

The method devised by Bauer et al. is able to agitate the phase space func- 
tion and build fluctuations, although their strength does not reproduce the 
expected value for fermionic systems. However, a careful modiflcation of the 
original procedure can considerably improve the results. 

The new procedure can be summarized in the following steps: 

(1) The choice of the two colliding partners closely traces the standard recipe. 
If the two test particles are allowed to collide, two clouds of Ntest particles 
will be moved, with conditions specified below. 

(2) Only particles within a sphere, in coordinate space, around the center of 
mass of the two partners ii and 12 can belong to the clouds. The distance 
criterion is |rj — TcmI^i, "^2)! < c^r, where dr is a free parameter (see later 
discussion). Then, for this considered space sphere, a grid is introduced 
in momentum space, around ii and ^2, the size of each cell being Vceii- 

(3) Given the momentum space cells / and J containing the partners ii and 
12, we consider the cells /' and J', corresponding to the final positions, 
in a rotated frame, as indicated in Fig. [2J For a given set of initial and 
final cells, the number of test particles that will be actually moved to 
final states is equal to the minimum between the occupation of the initial 
cells, n, and the availability of the final ones n' = {Nmax — n'): 

ni(/,/'; J, J') = mm{ni,nj,nr,nj>), 

where Nmax = gVrVceii/h^ is the maximum number of test particles that 
can stay inside a cell. It should be noticed that this choice corresponds 
to the maximum of possible moves. 

(4) Surrounding momentum space cells are searched with the same prescrip- 
tions, until two entire nucleons are found. The search procedure is sym- 
metric with respect to the center of momentum of the two partners and 
random in the vicinity of the original cells (I, J). 

(5) Finally, the two clouds are moved to the new states of the rotated frame 
(see Fig. 2), so that the conservation of energy and momentum is auto- 
matically obtained. However, a further check is performed. In fact, due 
to the finite number of test particles, the one-body density is not per- 
fectly homogeneous inside each cell. This causes slight violations of the 
conservations laws. Hence the origin of the momentum space grid is even- 
tually slightly displaced in order to have a perfect energy and momentum 
conservation. 

This method involves two parameters, namely the radius of the sphere in 
r space, dr, and the size of the momentum cells, Vceii- The radius dr fixes 
the spatial extension of the nucleon, that in turn infiuences its spreading in 
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Fig. 2. Schematic picture of the cohisional procedure in the improved pseudo-particle 
correlation method. 

momentum space. Hence this parameter fixes somehow the extension of the 
nucleon wave packet in phase-space and could affect the transport dynamics. 
It can be constrained by physical arguments, and in general it should depend 
on the physical properties of the system, such as its dilution. 

4.1 General comments 

The method proposed here has several advantages with respect to the original 
one of Ref. [23J. It is still very flexible, and can be easily implemented into 
existing transport codes; moreover, it is quite fast, and does not require huge 
amounts of resources. The Pauli blocking is locally checked in each collision, 
so in principle this procedure satisfies the basic requirements to preserve the 
average evolution and to yield the expected fluctuations. Besides, differently 
from the original method, where the nucleon is essentially spherical in phase 
space, here it is allowed to assume any shape in momentum space, since the 
only restriction is given by the cell size Vceii- In other words, it is possible 
to take into account momentum deformations of the nucleon wave packet, 
due to the momentum distribution of the system in particular conditions. For 
instance, when equilibrium is approached, the Pauli principle allows collisions 
essentially for nucleons lying on the surface of the Fermi sphere; therefore such 
nucleons should not have a spherical shape in momentum space. 

It is also interesting to make a connection between the proposed method and 
the lattice simulation method to build fluctuations. 

The numerical procedure described in Ref. [18], where the phase space is 
discretized by a fixed grid of cells, allows to construct the correct value of the 
fluctuation amplitude in each cell of the grid, generally taken of volume 
(but in principle it can be even smaller), on the basis of the fluctuating gain 
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and loss terms. 



On the other hand, our procedure is based on a different philosophy: Instead 
of constructing fluctuations in fixed cells, when a collision occurs we move 
all together the phase space volume containing one nucleon. The shape of 
this volume is not fixed a priori and depends on the location of the four 
points (pi, P2, Ps, P4) that have been chosen as initial and final centroids of 
the considered two-body collision. This means that the "correlation volume" , 
i.e. the volume of the sphere that envelops all test particles that move all 
together in each collision, is not fixed, but has to be considered as a dynamical 
variable that generally exceeds and could even reach rather large values. 
In other words, compared to what is done in Ref. [18], we build fluctuations 
in volumes containing one nucleon, but these volumes do not correspond to 
single cells of a fixed grid: instead, their contour changes from one collision to 
another, according to the prescription adopted to construct the clouds around 
the two first colliding test particles ii and 12- 

But how can the two proposed methods be connected to the original BL 
equation ? It should be noticed that the Boltzmann-Langevin equation ([T]) 
is a semi-classical equation, derived in the "continuum limit" (/i ^ 0). This 
implies that the evolution of the system under consideration is completely de- 
termined by a smooth distribution function, i.e. it behaves as a fluid in phase 
space. In other words, the volume As of the phase space cells for which the 
fluctuating transition rate is evaluated, is supposed to be much larger than . 
The Boltzmann-Langevin equation is essentially based on the idea of a fluc- 
tuating collision rate among phase space cells containing many nucleons. The 
formalism does not provide any additional information about the "structure" , 
i.e. the shape of the single nucleon wave packet. 

In both numerical implementations described above, two-body collisions among 
nucleons are treated as a stochastic process. Only the definition of nucleon 
wave packet and/or of the elementary cells where fluctuations are built change 
from one method to the other. As discussed above, this ingredient is not 
contained in the BL equation and such information is beyond its derivation. 
Hence, both procedures can be considered as correct implementations of the 
Boltzmann-Langevin theory in the "continuum limit", the definition of "nu- 
cleon wave packet" being rather arbitrary. 

However, it is also very interesting to investigate the differences between the 
various possible methods used to build the nucleon clouds, that in turn in- 
fluence the correlation volume. In particular, one can also try to force the 
procedure in order to reduce the wave packet smearing and to correlate par- 
ticles inside smaller volumes (close to /i^). 
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4-2 Details of the model 



In the present first implementation of tlie metliod, similarly to what was done 
for the lattice simulations of Ref.[T2], we neglect the evolution in coordinate 
space, extending dr to the whole coordinate space. The cell size Vceii is then 
constrained by the following arguments: it should be small enough to allow an 
accurate check of the Pauli blocking, but large enough to contain a sufficient 
number of test particles to reduce numerical uncertainties. In any case, Vceii has 
not to exceed the volume where, at most, one nucleon can be accommodated: 
Vp = h^/{gVr), being Vj. the volume in the coordinate space (K- = Ittc?^). In 
the following we take g = I, unless expressly indicated. Since we neglect the 
coordinate space, we will refer only to volumes V in momentum space. The 
expected equilibrium value of the variance aj in such volumes is: 

aj=/e./e.^ = ^ (14) 

where Ny is the number of nucleons that can be accommodated, at most, in the 
considered volume V. It should be noticed that, according to our procedure to 
build the nucleon wave packet, the expected value of fluctuations, (Eq JT^ . can 
be reproduced only in volumes V larger than the average correlation volume, 
while it is underestimated in smaller cells, since fluctuations are built on a 
larger scale. 



4-3 Illustrative results 



In order to illustrate how our procedure works, we will first consider the very 
simple case of particles moving only along one direction. We take a large 
system, containing 1000 nucleons with average occupancy -< f >~= 0.5. So 
the total extension in momentum space of the system is 2000 x Vp. Then 
the system is divided into = 4000 smaller cells that can contain at most 
half nucleon, whose occupancy can be either or 1. Hence Vceii = Vp/2. The 
occupancy is randomly chosen at the beginning and collisions are performed 
until equilibrium has been reached. The initial {pi,P2) and final (^3,^4) states 
of a collision are chosen randomly among the N cells. We neglect energy and 
momentum conservation in this simple example. Moreover, for simplicity, we 
take P2 = —pi and p^ = —p^. Hence only half of the considered space is 
independent and, in this situation, the probability for a collision to happen 
is proportional to the product /(pi)(l — fips)). According to Eq. f[T^ . the 
expected fiuctuation value, in a given momentum space volume V, will be 
(Tq = 0.25 Vp/V. The general prescription (see Sect. H]) states that when a 
collision happens, a full nucleon has to collide. Therefore, we move together 
the cell i and the closest occupied cell that finds an empty cell in the final 
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position. This example can be seen as a simple implementation of the general 
procedure when considering nucleons represented by only two test particles. 
One can easily realize that the "correlation radius" , i.e. the distance between 
the two cells that move together, changes from one event (collision) to another: 
sometimes the two co-moving cells are neighbours, sometimes they can be 
rather far. Of course, for practical reasons, we have to restrict our search within 
a given radius. For instance, considering 10 cells on each side (with respect to 
i), we are able to reconstruct the nucleon and to perform the collision in more 
than 98% of the cases. 
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Fig. 3. Variance of the distribution function as a function of the number Ny (see 
text), in the idealized case of particles moving only along one direction. 

In Fig. [3] we present the results for the fluctuation variance as a function of 
the number Ny of "nucleon" cells (of volume Vp) that are contained inside the 
volume V where the function / is evaluated. One thousand events have been 
considered. The fluctuation value presented in the Figure is multiplied by Ny. 
As one can see in the Figure, the correct value of the fluctuations, i.e. NvcTq = 
0.25, is approximately recovered (within ~ 5%) for volumes V ^ 20 Vp. The 
slight deviation is due to the upper bound considered for the "correlation 
radius" and to the constraints imposed by the total mass conservation, that 
are also responsible for the decrease of fluctuations observed for large Ny 
values. 

From this simple example, it already appears that the exact fluctuation value 
(iVyCJo = 0.25) is reached only asymptotically. This is due to the arbitrariness 
of the nucleon shape and of the centroid position. Hence, when testing a given 
procedure to solve the stochastic coHision integral, the resulting fluctuations 
have to be compared to the expected values in large volumes, in which all 
different possibilities (nucleon configurations) may be accommodated. 

These considerations would apply also to the original procedure proposed by 
Bauer et al. Thus, in Ref.[29], the comparison between the obtained fluctua- 
tion variances and the expected ones should have been done considering also 
volumes containing many nucleons, and not only cells of volume Vp where, 
as seen in Figl31 fluctuations are naturally suppressed by the fact that the 
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nucleoli wave packet may take any shape. However, this does not affect the 
conclusions drawn in [29], mainly based on the wrong profile of the density 
variance (and not on its absolute value). 

In the next Section we will move to describe more complex geometries and 
different procedures to build the nucleon wave packets. 



5 Fluctuations of a Fermi gas: Analysis on the Fermi surface 

Here we will discuss our procedure to build fluctuations in the context of a 
system of fermions having a given density and temperature, i.e. initialized ac- 
cording to a Fermi-Dirac distribution. This is a situation that is easily reached 
in the course of a dissipative nuclear reaction, after the initial coUisional shock 
[8]. Since collisions happen mostly on the Fermi surface, let us first consider 
only particles within a thin stripe dp around the Fermi momentum, where 
^/^=0.5. 

For those colliding particles the relative velocity Vi2 is constant, since only 
particles with opposite momenta can collide, in order to stay on the Fermi 
surface. So the transition rate is essentially governed by the product /1/2/3/4. 

This simplified system will allow us to study the sensitivity of the results to the 
ingredients of the procedure used to build, step by step, the nucleon clouds. 
As also discussed before, we will show that the correct fluctuation amplitude 
is recovered in the "continuum limit", i.e. in large volumes. However, in the 
following we will mostly concentrate on the fluctuation amplitude in volumes of 
size Vp, in order to probe the extension of the nucleon wave packet. The latter 
could influence significantly the evolution of actual nuclear collisions, where it 
could be important to work with more compact nucleon configurations. 

Using this simple model we will find that building the maximum fluctuation 
value (e.g. cTq = 0.25) in such a volume, i.e. correlating particles inside Vp, is 
a quite difficult task and additional efforts are required in order to reduce the 
nucleon cloud smearing. 

5. 1 Details of the simplified model 

In order to keep the particles on the Fermi surface, the original procedure 
described in Sect. H] was slightly modified so that the modulus of the mo- 
mentum of the particles does not change: only particles belonging to cells 
with opposite momenta can collide. Then, for symmetry reasons it will result 
/(Pi) = /(P2 = ~Pi)- The transition probability finally depends only on two 
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values of the distribution function, /(pi) and /(ps). We use spherical coordi- 
nates (^cos{6) , (f)j as independent variables. The collision partners of particles 

in the cell cos{6), (f)j will be the ones in the cell cos{it — 6), iT+(f)j . The choice 
of these coordinates implies the grids used in the procedure to be fixed along 
the cos(^) axis. They can slide along the axis, where periodic conditions are 
taken into account. 

The procedure is as follows: a test particle ii is chosen at random (together 
with the corresponding partner ^2) and a grid is constructed, centered around 
ii and ^2 (see Fig. 2). The rotation angles, also taken from a flat distribution, 
fix the final states and the new frame corresponding to these final positions. 

We consider a shell with 40 x 40 cells, whose size is chosen to accommodate 
one nucleon at most, so that Vceii = Vp. We have checked that the results do 
not depend on the number of cells, i.e. the number of nucleons considered, 
by performing the same calculations on a 64 x 64 system. We use 500 test 
particles per nucleon. 

The growth of fluctuations is described by the variance aj of the distribution 
function, as discussed before. Usually it is evaluated performing an ensemble 
average over a large number of events having the same initialization. However, 
in this scheme the mean value -< f >- does not depend on the position p and 
is not changing in time. Considering the large number of cells employed, we 
can also follow the evolution of the system by calculating in a single event the 
following quantity: 

^/ = ((/^-(/))'U- (15) 
where (■) denotes the average with respect to the cells and (/) = 0.5. 

5.2 "Fixed grid" calculations 

In order to understand how the recognition of the nucleon wave packet can 
be influenced by the details of our procedure, we first further simplify the 
problem by using a single fixed grid to make the collisions, i.e. test particles 
are placed at (and can only move to) the center of the cells of the grid. In this 
case we have a single reference frame, so the effect of a collision (in the space 
(cos(e),0) ) is a mere translation to other states of the same grid. Apart 
from the use of test particles, and related differences in the method followed 
to solve the collision integral, in this way the procedure becomes similar to 
the one used in Ref. [T8] . 

We will take, as initial conditions, / strictly equal to 0.5 in each cell of the 
fixed grid, corresponding to 800 nucleons. By construction, once collisions are 
performed, the only possible values for / are 0, 0.5 and 1. 
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The maximum fluctuations in Vp {a^ = 0.25) would be reached only if, in each 
cell, / is either zero or one. Achieving this result is not a trivial task: Cells 
starting with / = 1 can be partially emptied, and cells with / = partially 
filled. 

The dot-dashed curve in FigjH upper panel, shows the variance of the distribu- 
tion function calculated in the cells of the fixed grid, with increasing number 
of collisions ( function of time). It saturates approximately when al- 

most all nucleons have suffered one collision. In the equilibrium configuration 
the number of cells with the three possible values of / is roughly the same, as 
shown in Fig. H] (lower panel), so the fluctuation variance is o"j ~ 0.175, which 
is about 2/3 of the maximum. 
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Fig. 4. Upper panel: variance of the distribution function in Vp, as a function of 
the number of collisions, for three different cases: fixed grid with f{t = 0) = 0.5 
(dot-dashed) , fixed grid with randomly distributed test particles (full) , moving grid 
(dashed). Lower panel: Distribution of the number of cells as a function of their 
occupation / at equilibrium, for the three cases indicated above. 
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The fluctuations reached in Vp also depend on the variance of the initial con- 
figuration. If we initialize the system distributing the test particles randomly, 
the occupation in the cells has an average value (/) = 0.5 and a width given 
by af{t = 0)/(/) = iV-^/^ where N is the average number of test particles 
per cell {N = 250). The distribution function can assume any value between 
and 1, as shown by the full curve in Fig. HI lower panel. The initial continu- 
ous distribution of values of / causes a reduction of the fluctuations in Vp (full 
curve in Fig. (jl]), upper panel) with respect to the previous case. As expected, 
the fact that / varies continuously between and 1 contributes to enhance 
the smearing of the nucleon wave packet and hence to reduce fluctuations in 
Vp. 

5.3 The centroid degree of freedom 

Now we turn to a more general situation, relaxing the condition of the fixed 
grid. Hence here we would like to test the full procedure, described in Section 
4, that we will finally follow in the full 3D case. The two reference frames 
used to make the collision (the one concerning the initial states and the one 
for the final states) are now shifted along the axis, since the angle of the 
initial and final states can assume any value between and 2n, according to 
the initial and final positions of the colliding test particles, ii and 22- 

This carries two important consequences. First of all, even being the Pauli 
blocking checked for each collision in a particular frame, the prescription / < 1 
can be violated in another frame. However, we checked that this violation 
concerns only about 5% of the cells. Second, the correlation between the cells 
representing the reconstructed nucleon, which is maximum in the search frame, 
will be systematically broken in another (shifted) reference frame, increasing 
the nucleon delocalization and thus reducing fluctuations in Vp. 

The dashed curve in Fig. H] shows the fluctuations obtained in Vp in this 
case. Since this feature is due to the arbitrary position of the centroid of 
the constructed nucleon with respect to a fixed reference frame, we will name 
it centroid effect. We notice a fluctuation reduction of ~ 60% with respect to 
the fixed grid case (dash-dotted line in the same Figure). 

5.4 Sensitivity to the dimension of the grid cell 

Up to now, the step of the grid has been chosen so that Vceii = Vp. However, in 
general, the size of the search cell has to be smaller than Vp, especially when 
the nucleon wave packet is rather delocalized in p space (being more compact 
in coordinate space). In fact, small search cells are required for the Pauli 
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blocking to be accurate. Hence we want to study the sensitivity of the results 
to the size of the search cell. Let us start with the fixed grid scheme. In Fig. [5] 
we plot the fluctuation variance obtained dividing the step in the direction, 
0step, by two or four; the volume of the search cell is then Vceii = Vp/2, 
respectively. This is equivalent to shift, for each collision, the origin of the 
grid by 0.5 or 0.25 (pstep respectively, and acts also as a sort of centroid effect, 
thus reducing the variance in V„ (see solid and dashed lines in Fig. [5]). The 
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Fig. 5. Variance of the distribution function in the fixed grid scheme and in the 
general case, for two sizes of the search cell. The lines on the right indicate the 
values for Vceii = Vp, also shown in Fig. [H 

value obtained with Vceii = Vp/2 is lower than the corresponding case with 
Vceii = Vp, indicated as "fixed grid" in the Figure. Moreover, it further reduces 
for Vceii = Vp/A. Nevertheless, this smearing effect is partly compensated by 
the fact that the procedure used to select the cells belonging to the nucleon 
clouds favours compact, symmetric configurations. 

It is interesting to evaluate the effect of the use of a small search cell on 
the fluctuations also in the general case, i.e. when the fixed grid condition is 
relaxed. Intuitively, we expect a small effect on the variance, as long as Vceii 
is not too small: the main reduction of fluctuations comes already from the 
centroid effect (see previous Subsection), which is equivalent to an infinitely 
small step. In fact, the differences between calculations with Vceii = Vp/2 and 
Vceii = Vp/A (Fig. [5l dotted and dash-dotted curves respectively) are smaher 
than in the previous calculation with fixed grid. The reduction with respect 
to the case with Vceii = Vp (denoted as "centroid effect" in the Figure) is also 
smaller. 

5.5 Conclusions on the search procedure 

In conclusion, it appears from our calculations that the smearing of the nucleon 
wave packet is rather dependent on the procedure adopted to build nucleons 
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and is particularly affected by the so-called centroid effect, related to the 
use of randomly distributed test particles for the phase-space mapping. As a 
consequence, fluctuations calculated in volumes Vp appear reduced. As shown 
above, the use of a flxed grid to locate the initial and flnal cells of the colliding 
test particles would help to increase the value of the variance in Vp. However, 
energy and momentum conservation would not be exact, due to the phase- 
space discretization. On the other hand, the standard test particle method 
allows a good mapping of the phase space and a good resolution of the average 
collision integral, with perfect energy and momentum conservation. However, 
when constructing fluctuations, the centroid effect may destroy the correlations 
already built and increase the nucleon cloud smearing. However, the expected 
statistical value of the fluctuation variance is recovered in large volumes. 

This last point has been checked for all the simplifled situations considered 
here (Fermi surface). For instance, this analysis is shown in Figl6]in the case of 
the flxed grid with Vceii = Vp and test particles randomly distributed inside the 
system (last case of Section 5.2). Here fluctuations are evaluated in volumes 
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Fig. 6. Variance of the distribution function (rescaled by A^v') as a function of the 
number Ny of nucleon cells Vp, obtained in the fixed grid case with Vceii = Vp and 
test particles randomly distributed. 



V, by considering an ensemble of 800 events. As expected, the fluctuation 
variance, rescaled by A'^v') almost approaches the expected value as the volume 
V increases. As explained above for the idealized case presented in Fig.3, 
the maximum value, Nya^ = 0.25, is not exactly reached because of the 
upper bound imposed to the correlation radius in the nucleon reconstruction 
procedure and of mass conservation constraints. 
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6 Non-equilibrium situations 



In the Fermi energy domain, two colliding nuclei may be schematically repre- 
sented by touching Fermi spheres in p space, so it is interesting to study the 
development of fluctuations starting from this initial distribution. 

In our simplifled 2D scheme, we can mimic the two spheres using a "chess- 
board" initialization, that is subdividing our momentum space in four re- 
gions where / is alternatively zero or one. Although the system is in a non- 
equilibrium situation, the variance, defined as in Eq. ffTSj) . starts from its max- 
imum value (since (/) = 0.5). It can be shown that in the fixed grid case the 
equilibrium configuration consists of a random distribution of full and empty 
cells, and the variance keeps constant in time. However, a very different result 
comes out either changing the search step or removing the fixed grid con- 
straint. In both simulations the variance rapidly decreases because the initial 
correlation is destroyed by the centroid effect. 

As an illustrative example, in Fig. [7] we plot aj for the choice Vceii = Vp/2, 
fixed grid. Comparison with the corresponding curve for the configuration 
/(p) = 0.5 ± af(t = 0) (dashed, see Section 5.4) shows that the two calcu- 
lations converge to the same value. We get a confirmation of these findings 
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Fig. 7. Comparison between variances obtained with the "chess-board" initialization 
(solid line) and with the case / = 0.5 ± af{t = 0) (dashed curve) in the fixed grid 
case for Vceii = Vpf^- 

by removing the fixed grid constraint. Also in this case, at the end we find 
no difference with the case described before (initialization with randomly dis- 
tributed test particles. Section 5.3), see Fig. [HI Moreover, in this case the 
convergence between the two calculations is rather quick. Then we can con- 
clude that the final outcome does not depend on the initial conditions of the 
system, but only on the final equilibrium situation, as it should be. 
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Fig. 8. Comparison between variances obtained with the same initializations of Fig. 
[7] but for Vceii = Vp and without the fixed grid constraint. 

7 Results for full 3D simulations 

Guided by the simplified schematic cases illustrated above, let us now discuss 
the full case of a nuclear Fermi gas of particles at a given temperature and 
density interacting through hard two-body scattering. 

As in the previous examples, since we focus on fluctuations in momentum 
space, only one large cubic cell is present in coordinate space, and all particles 
can be chosen to collide (no restrictions in r space). The size of the box is 
L = 26 fm, and we consider 2820 nucleons, so that the density has the satu- 
ration value Qo = 0.16 fm~^; each nucleon is represented by a collection of 500 
test particles. Besides, we do not consider any distinction between neutrons 
and protons, so that one nucleon occupies at least a phase-space volume h^/A 
((7 = 4). We initialize the momenta so as to reproduce a Fermi-Dirac proflle 
corresponding to a temperature of 5MeV. Finally, we consider a constant 
cross section a = 160 mb. In these calculations the volume Vceii corresponds to 
a cube of side Ig = 30 MeV/c (and coincides with Vp = h?/{AL'^)). We notice 
that, for the full Fermi g grid in Cartesian coordinates is easier to 

use, with respect to spherical coordinates, in the nucleon construction proce- 
dure. However, spherical coordinates will be adopted to analyze the resulting 
fluctuations. Calculations are stopped when the fluctuation variance saturates. 

First of all, we checked the effect of our fluctuating collision integral on the 
average evolution of the system. In Fig. we plot the energy proflle of the 
distribution function at the initial time (solid) and at the time when the 
fluctuation variance saturates and we stop the calculation (dashed line) . Slight 
changes are observed, probably due to the flnite extension of the nucleon 
packet, that induces a discretization of the phase space. 
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Fig. 9. Distribution function profile at tlie initial time t = (solid) and at the time 
t = fOO/m/c (dashed line). 

We also compared the average collision rate with analytical expectations: we 
find that the number of collisions is slightly larger than the analytical one, 
but differences are within 10%; on the other hand, similar deviations are also 
found in usual BUU implementations, due to the approximate mapping of 
phase space induced by the finite number of test particles. 

In the following we analyze in more detail the capability of the procedure 
to build the fluctuations in the considered 3D scenario. We investigated the 
region in p space where collisions are most probable and the extension of the 
nucleon cloud. From the left panel of Fig. (TU] it is evident that most collisions 
involve nucleons lying on the surface of the Fermi sphere, as expected. Besides, 
the right panel illustrates the extension of the clouds in the radial direction, 
calculated according to the definition: 

Ap = ^I^Jjh^'WU^d, (16) 

where (p) denotes the center of momentum of the cloud and Pi is the modulus 
of the momentum of the test particles contained in the cloud. 

From this distribution one can deduce that the clouds correspond to a rela- 
tively narrow stripe, 2Ap 35MeV/c, in the region where / is different from 
either zero or one. Moreover, we have checked that the extension of the con- 
structed wave packet is similar in the three directions Px,Py,Pz, as expected. 
We notice also that the tail of the distribution (Fig. 10, right) extends up to 
large Ap values, ~ 50MeV/c. Although the nucleon clouds extends, on av- 
erage, over a volume {2Ap)^ ~ {35 MeV/ c)^, that is not much larger than 
Vp, the centroid effect, i.e. the fact that the nucleon centroid may be located 
anywhere, causes a significant smearing of fluctuations, as already discussed 
in Section 5. In fact, we find that, for our Fermi gas calculations, fluctuations 
in Vp are considerably underestimated, as shown in FigJTTl where the variance 
(j| is reported as a function of the energy E = p'^/{2m). In the following, 
to test our procedure, we will evaluate fluctuations considering big volumes. 
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Fig. 10. Left panel: distribution of the position of the centroid of the nucleon 
"clouds" in the \p\ direction. Right panel: distribution of the size of the nucleons on 
the same axis. 
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Fig. 11. Profile of the fluctuation variance, as a function of the energy E = / (2m), 
in cubic cells of size l^, for a Fermi gas at equilibrium. 

that contain many nucleons. However, it should be noticed that, in the geom- 
etry considered, it is not trivial to predict which is the size and shape of the 
most appropriate volume to recover the expected value of fluctuations. While, 
from one side, it is desirable to consider volumes containing a large number 
of nucleons, from the other side one should keep in mind that spurious fluc- 
tuation reduction, due to the conservation of the total number of particles, 
may be a problem. Hence we have performed a systematic analysis of fluctua- 
tions as a function of the coordinates (p, 6), integrating over the angle 0. The 
corresponding volume employed for the calculation of the variance is given by: 

= 27r^Acos^ 

We consider volumes corresponding to flxed angular spreads 9 and/or to a 
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given step Ap^. In Fig. [12] we fix the angular spread, 6 = 30°, and we con- 
sider four different steps, Ap^, corresponding to different volumes, i.e. to a 
different number of nucleons Ny that can be contained inside. The consid- 




Fig. 12. Fluctuation variance as a function of the energy, for a fixed angular spread 
e = 30° and different steps Apstep = (Ap3)i/3, 

ered volumes can contain up to a rather large number of nucleons (about 
90). The variance, multiplied by Ny, is plotted as a function of the energy 
E. We find that the peak of the variance is around the Fermi energy, as ex- 
pected. The rescaled variance gets larger when we increase the step Ap^. How- 
ever, even in the best situation, the expected equilibrium value at the peak 
{NyaQ = 0.25) is underestimated by a factor ~ 2.3. Therefore, for a fixed 
value of Apstep = {Ap^Y^^ = 190 MeV/c, we tested different angular spreads. 
In Fig|T3]we present the obtained results, rescaled by Ny. Now we can observe 
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Fig. 13. Fluctuation variance as a function of the energy, for a fixed 
Apstep = 190 MeV/ c and different angular spreads (curves with circles). The thick 
full and dashed lines are explained in the text. 

that, increasing the angular spread, the variance first increases, then decreases 
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again, due to the mass conservation constraints. The variance gets larger when 
considering the angular spread 9 = 20°, but the expected value is still under- 
estimated by about a factor 2, at the peak. This result can be understood 
on the basis of the results obtained in Vp and displayed in FigJTTl In fact, in 
order to recover the expected fluctuation value, the volume employed for the 
calculation of the variance must be much larger than V^^^ in all the directions 
{p,cos{6),(j)). Due to the considered geometry, while this can be achieved in 
the angular coordinates, this cannot be easily obtained along the p direction 
without losing locality for the value of / and related fluctuations. Actually, 
in the case of Ap^tep = 190 MeV/c, the corresponding extension along the p 
direction, Apmod, of the volume considered to evaluate the function / is ap- 
proximately equal to 30MeV/c = Vp^^ at p = pp- We thus expect fluctuations 
to be reduced, with respect to fegfeg, by a factor a^^^ ^ 0.48, being a ~ 0.11 
the reduction factor in Vp, as evaluated from Fig{TTl The expected fluctuation 
reduction can be extracted also from the schematic calculations described in 
Sect. 5. In fact, from the inspection of FigSJ it appears that, in the flxed grid 
scheme, fluctuations in Vp are reduced by a factor 0.56. Relaxing the flxed grid 
constraint in one dimension, fluctuations are further reduced by a factor 0.64 
with respect to the flxed grid case. We thus expect a fluctuation reduction 
by a factor 0.56^/^ x 0.64 ^ 0.5 (at the Fermi surface). This is close to our 
result at E = Ep = p\/{2'm) (see the dashed line with circles in Fig. 13, for 
instance). 

One may try to introduce this correction in the results of Fig. 13. In the 
considered geometry, the reduction factor a of the fluctuations in Vp depends 
on the energy E, i.e. one can write, for the variance in Vp-. 

aj{E,Nv = l) = E{E)a{E), (17) 

where F(E) indicates the expected fluctuation value, that in our case should 
coincide with fegfeg- 

By approximating the volume dependence of the rescaled variance by a linear 
behaviour for small volumes (see FigOfor instance), the fluctuation variance 
presented in Figs. 12-13 can be rewritten as: 

NyajiE, Ny) ^ F{E) a"\E)^. (18) 

Vp 

We notice that Ap^od is also depending on E. 

Eqs.f lT7|)ffTgl) allow to extract the suppression factor a{E), as well as the func- 
tion F{E). The latter is displayed in Fig. 13 (thick dashed line). One can see 
that the expected fluctuation value, fegfeg (thick solid line), is well repro- 
duced at the Fermi surface, while it is underestimated, within 30%, expecially 
for outer regions. This underestimation indicates that in our procedure, due 



26 



to the finite extension of tlie nucleon cloud, tlie regions far from the Fermi 
surface are less involved in the building of fluctuations than they should be 
and can be considered as an intrinsic limitation of the method. 

However, the overall shape of the fluctuation variance is reasonably repro- 
duced, relative to //: This is a nice indication that the Pauli blocking is 
effective, and the equilibrium distribution is almost unchanged by the proce- 
dure. Therefore, it really represents a remarkable improvement with respect 
to the original method by Bauer et al. (see Ref. [29]), which makes us confident 
about applications to nuclear reactions. 

We can directly visualize the effect of the fluctuating term by selecting a thin 
region in momentum space around the Fermi momentum. For this purpose, 
we adopt a new set of coordinates {dp, pd9, psin(6')(i0); for fixed p and dp, the 
distribution function depends only on two coordinates, namely: 

f(p,e,^)^ fie,r]) 

with 7] = sm{6)(j). This set of coordinates correspond to displacements of 
equal length along the three directions: p,0,(f). We choose p = 2Q0 MeV/c 
(approximately equal to the Fermi momentum) and dp = 10 MeV/c, and in 
Fig. [m we plot the distribution function f{6,i]) at two different times (initial 
and final times). We recognize the sinusoidal profile on the rj axis due to the 
modulation given by the sin(^) factor. At the initial time, / is nearly uniform, 
and its fluctuations are simply due to the numerical noise induced by the finite 
number of test particles. At later times, we observe the growth of fluctuations, 
evidenced by the typical structure with "peaks and holes" . 




Fig. 14. Distribution function f{9,r]) at time t = (left) and t = 100 fm/c (riglit). 
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8 Methods to increase fluctuations in Vp 

In actual nuclear collisions the extension of the nucleon wave packet in phase- 
space may affect the transport dynamics and it may be important to correlate 
the test particles belonging to the nucleon clouds inside smaller volumes. In 
the context of the simplified model (Fermi surface) discussed in Sect. 5 we 
attempted to find a method to reduce the nucleon cloud smearing and enhance 
the fluctuations in Vp without losing the consistency of the original procedure. 

First of all, it should be remembered that the collision probability is given by 
/(Pi)/(P3) (see Sect. 15.11) : therefore, we need to consider only two locations 
in p space. The main idea is to force the system to fully empty the cells 
associated with the initial state or fully fill the final state when a collision 
occurs. This condition is not automatically achieved for all cells contributing 
to a given collision. In fact, the number of particles, nt{i), taken from a cell i 
and employed to build the nucleon obey to the constraint: 

= Ntest 

i 

being Ntest the number of test particles per nucleon. Hence the number of 
test particles that are really taken from a cell is not always equal to rit = 
'min{ni,n3), but instead to n[ = min{nt,nres), being 

i-l 

i'=l 

The idea is to try to avoid the situation when the number of particles moved is 
less than n^. In this way fluctuations get to a larger value, since we favour the 
transitions where cells are completely filled or emptied. In order to implement 
this idea, when constructing the nucleon clouds we choose, among the possible 
ones, the cells for which the probability = ^ is large. 

The outcome of this method is plotted in Fig. [TB] for the fixed grid case, where 
the result of Fig. H] is also shown for comparison (solid line). As the dashed 
curve in Fig. [15] clearly shows, this optimization procedure is successful in 
enhancing the fluctuations, which eventually reach almost their upper limit 
o"o = 0.25. Actually, this maximum value is not obtained due to the fact 
that only adjacent cells are considered in the nucleon construction. In this 
way, some half-filled cells remain "isolated" and are not involved in a collision 
anymore. Indeed, at the time when we stop the calculation, isolated cells 
amount to about 5% of the total. 

The method works also in more general situations, for instance in the case 
corresponding to fluctuating initial conditions (see Sect. 5. 2). The comparison 
between variances with and without the optimization procedure is shown in 
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Fig. 15. Variance of the distribution function in Vp, as a function of the number of 
collisions, in the fixed grid scheme. Solid line: no optimization; dashed line: opti- 
mized procedure. 




Fig. 16. Variance of the distribution function in Vp, as a function of the number of 
collisions, in the fixed grid scheme with fluctuating initial conditions. Solid line: no 
optimization; dashed line: optimized procedure. 



Fig. [161 As in the previous case, the optimization favours a faster growth 
of the fluctuations. Their saturation value for y = Vp is now enhanced, by 
~ 20% with respect to the standard procedure. This is an indication that the 
procedure is able to increase the number of empty and full cells. This feature 
is also evident from the distribution of values of / plotted, for the same two 
simulations of Fig. 16, in Fig. [T71 

In the end, this procedure helps the nucleon packets to keep a similar, more 
compact shape in p space during the time evolution of the system, from one 
coUisional event to another. As a consequence, nucleons are more localized 
and fluctuations in Vp appear enhanced. 
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Fig. 17. Distribution of the number of cells as a function of their occupation / for 
the same two calculations of Fig. [TBI 

9 Conclusions 

We have proposed a new numerical implementation of the full Boltzmann- 
Langevin equation in 3D. The stochastic character of the two-body collision 
integral, that in standard transport codes applies only to the single test parti- 
cles, leading to a strong reduction of fluctuations, is recovered, at the nucleon 
level, by a careful reconstruction of the nucleon wave packet once collisions 
occur. This is achieved by moving, in each collisional event, two entire test 
particle clouds, corresponding to one nucleon each, as proposed in Ref . [23] . 
However, while in the original procedure the Pauli blocking was checked only 
for the centroids of the nucleon clouds, here we propose a new method that al- 
lows to construct nucleon wave packets by carefully checking the Pauli blocking 
for the entire swarm of test particles. The last point is essential in the con- 
struction of the correct fluctuation value in fermionic systems. The nucleon 
wave packet may take, in principle, any shape. Thus the correlation volume, 
i.e. the volume of the sphere that contains all test particles that move together 
is generally larger than . Moreover, the nucleon centroid is chosen randomly 
among the test particle distribution, leading to a partial overlap of nucleon 
configurations from one collision to another. All these effects lead to a reduc- 
tion of fluctuations in the cells of volume h^. However, the correct fluctuation 
value is recovered in larger volumes, where all possible nucleon configurations 
may be accommodated. 

The procedure, and its numerical ingredients, has been tested carefully in sev- 
eral idealized situations, studying in particular the relation between the results 
obtained for the fluctuations in cells of volume and the ingredients of the 
model. Considering a fermionic system at equilibrium, at a given density and 
temperature, we find that our procedure builds, within a good approxima- 
tion, the expected profile of the fluctuation variance, feqfeq, as a function of 
the nucleon energy. This can be considered as a significant improvement with 
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respect to the original procedure [23], that opens interesting possibihties of 
apphcations to nuclear collisions. In fact, the method proposed can be easily 
implemented into existing transport codes, and could allow to treat, in a more 
complete way, phenomena where large fluctuations or bifurcations occur and 
fluctuations are important. 

Finally we stress that the results presented here are relevant not only for 
nuclear fragmentation studies, but in general for the dynamical description of 
quantum many-body systems. 
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